Selection of Ideal Reference Genes for Gene Expression Analysis in COVID-19 and Mucormycosis

ABSTRACT Selection of reference genes during real-time quantitative PCR (qRT-PCR) is critical to determine accurate and reliable mRNA expression. Nonetheless, not a single study has investigated the expression stability of candidate reference genes to determine their suitability as internal controls in SARS-CoV-2 infection or COVID-19-associated mucormycosis (CAM). Using qRT-PCR, we determined expression stability of the nine most commonly used housekeeping genes, namely, TATA-box binding protein (TBP), cyclophilin (CypA), β-2-microglobulin (B2M), 18S rRNA (18S), peroxisome proliferator-activated receptor gamma (PPARG) coactivator 1 alpha (PGC-1α), glucuronidase beta (GUSB), hypoxanthine phosphoribosyltransferase 1 (HPRT-1), β-ACTIN, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) in patients with COVID-19 of various severities (asymptomatic, mild, moderate, and severe) and those with CAM. We used statistical algorithms (delta-CT [threshold cycle], NormFinder, BestKeeper, GeNorm, and RefFinder) to select the most appropriate reference gene and observed that clinical severity profoundly influences expression stability of reference genes. CypA demonstrated the most consistent expression irrespective of disease severity and emerged as the most suitable reference gene in COVID-19 and CAM. Incidentally, GAPDH, the most commonly used reference gene, showed the maximum variations in expression and emerged as the least suitable. Next, we determined expression of nuclear factor erythroid 2-related factor 2 (NRF2), interleukin-6 (IL-6), and IL-15 using CypA and GAPDH as internal controls and show that CypA-normalized expression matches well with the RNA sequencing-based expression of these genes. Further, IL-6 expression correlated well with the plasma levels of IL-6 and C-reactive protein, a marker of inflammation. In conclusion, GAPDH emerged as the least suitable and CypA as the most suitable reference gene in COVID-19 and CAM. The results highlight the expression variability of housekeeping genes due to disease severity and provide a strong rationale for identification of appropriate reference genes in other chronic conditions as well. IMPORTANCE Gene expression studies are critical to develop new diagnostics, therapeutics, and prognostic modalities. However, accurate determination of expression requires data normalization with a reference gene, whose expression does not vary across different disease stages. Misidentification of a reference gene can produce inaccurate results. Unfortunately, despite the global impact of COVID-19 and an urgent unmet need for better treatment, not a single study has investigated the expression stability of housekeeping genes across the disease spectrum to determine their suitability as internal controls. Our study identifies CypA and then TBP as the two most suitable reference genes for COVID-19 and CAM. Further, GAPDH, the most commonly used reference gene in COVID-19 studies, turned out to be the least suitable. This work fills an important gap in the field and promises to facilitate determination of an accurate expression of genes to catalyze development of novel molecular diagnostics and therapeutics for improved patient care.

can be affected by several factors not only due to SARS-CoV-2 infections. For example, the age and health conditions of the subjects in the groups can be different. 3. Authors did not mention the quantity of cDNA was used in the RT-PCR reactions. 4. It was mentioned that the groups were categorized by institutional definitions, but the categorization of the groups must be approved by certain international bodies otherwise, it is hard to apply the findings of this study in general. 5. In Fig.2, some of the ct values for the GUSB or PGC-1a or GAPDH in pre-CAM or post-CAM are above ct 30.0, in some cases even above ct 34.0. I think the ct values above 30 is not reliable for calculations. Moreover, the PCR was performed at annealing temperature of 60 degrees for all the reference genes. However, the primers could have better amplification efficiency at other temperatures as well. Not sure, if authors have performed gradient PCR and selected the annealing temperature for the primers. 6. A 260/280 ratio of ~2.0 is generally accepted as "pure" for RNA. In Table 1, the RNA quality for the groups of Severe, pre-CAM and post-CAM is below 2.0 and varies a lot compared to other groups. This is very important factor in the whole experiment as the downstream steps depends on the purity of RNA used for cDNA preparation. 7. To support the RT-qPCR data, authors should have analyzed the expression of these reference genes by Western blotting from the PBMC's collected. 8. Based on normalization with CypA, NRF2 expression level did not correlate with the disease severity in the groups. Therefore, authors should have analyzed for at least 4-5 genes that may have differential expression in SARS-CoV-2 infections to prove that CypA is the suitable reference gene compared to GAPDH.
Reviewer #2 (Comments for the Author): The manuscript by Kumar et al., entitled "Selection of ideal reference genes for gene expression analysis in COVID-19, and Mucormycosis infection" identified among 9 different housekeeping genes, the CypA gene as the most suitable reference gene, which accurately captured the heterogeneity of infection and yielded a stable expression across different conditions of COVID-19 severity and Mucormycosis infection (Covid-associated mucormycosis -CAM). Briefly, the authors tested the expression stability of nine different housekeeping genes (including TBP, CypA, B2M, 18S, PGC-1α, GUSB, HPRT-1, β-ACTIN and GAPDH) in PBMC isolated from patients with different COVID-19 severity and CAM. The analysis of housekeeping genes during SARS-CoV-2 spectrum disease and CAM infection showed significant expression variations in most of candidate reference genes intra and inter-group. Only TBP, CypA, B2M and HPRT-1 had a comparable mean Ct value across SARS-CoV-2 spectrum and CAM infection. Then, the authors employed statistical algorithms to determine the stability indexes of each reference gene. GHPDH and PGC-1α had more deviations. Nonetheless, CypA is identifies as the most suitable reference gene in and CAM infections. Strength Studying the co-infection of SARS-CoV-2 and mucormycosis is relevant based on existing public health issue and the lack of studies published in the field. Also, the selection of accurate reference genes in different disease outcomes could advance the field.

Weaknesses
The observations postulated in the paper are novel and interesting. However, the discussion is not robust enough to make convincing conclusion. The authors should show more data from published studies related to their results.
Below are comments to strengthen the manuscript.  Figure 2 shows difference between healthy and Covid19 subjects with moderate or severe or CAM infections. GAPDH showed variations inter-groups, in which expression was comparable in health, mild and asymptomatic patients, while moderate and asymptomatic as well as Pre/Pos CAM had other expression patterns. How could you explain this pattern of high GAPDH variations inter-group ( Figure 2) whether in Figure 3b the data shows GAPDH as a reference gene more stable than CypA? 2) Figure 4: As the health subjects represent the baseline NRF2 expression, is there significative difference in NRF2 expression among samples from moderate, asymptomatic and/or pre/post-CAM diseases when CypA was used as internal control gene? If so, how can the NRF2 expression variability among them using CypA as internal control gene be explained?
3) Also, at line 164, the authors comment that NRF2 is suppressed in Covid-19 patients, which contradicts the Figure 4 data that shows NRF2 expression seems to be higher mainly in Moderate COVID-19 and pre-CAM disease comparing to health subjects even using the CypA as the internal control. Thus, the authors should specify the difference in the expression profile of NRF2 inter-group based on the usage of CypA.

Preparing Revision Guidelines
To submit your modified manuscript, log onto the eJP submission site at https://spectrum.msubmit.net/cgi-bin/main.plex. Go to Author Tasks and click the appropriate manuscript title to begin the revision process. The information that you entered when you first submitted the paper will be displayed. Please update the information as necessary. Here are a few examples of required updates that authors must address: • Point-by-point responses to the issues raised by the reviewers in a file named "Response to Reviewers," NOT IN YOUR COVER LETTER. • Upload a compare copy of the manuscript (without figures) as a "Marked-Up Manuscript" file. • Each figure must be uploaded as a separate file, and any multipanel figures must be assembled into one file. For complete guidelines on revision requirements, please see the journal Submission and Review Process requirements at https://journals.asm.org/journal/Spectrum/submission-review-process. Submissions of a paper that does not conform to Microbiology Spectrum guidelines will delay acceptance of your manuscript. " Please return the manuscript within 60 days; if you cannot complete the modification within this time period, please contact me. If you do not wish to modify the manuscript and prefer to submit it to another journal, please notify me of your decision immediately so that the manuscript may be formally withdrawn from consideration by Microbiology Spectrum.
If your manuscript is accepted for publication, you will be contacted separately about payment when the proofs are issued; please follow the instructions in that e-mail. Arrangements for payment must be made before your article is published. For a complete list of Publication Fees, including supplemental material costs, please visit our website.
Corresponding authors may join or renew ASM membership to obtain discounts on publication fees. Need to upgrade your membership level? Please contact Customer Service at Service@asmusa.org.
Thank you for submitting your paper to Microbiology Spectrum.

Comments
The manuscript by Sunil et al., has carried out very important study during this ongoing COVID-19 pandemic. Although, this study provides strong rationale for identification of appropriate reference genes in the viral infections but has some serious problems in the experiments and result interpretations. Overall, the results are not sufficient and convincing to their claims. Here are my comments for the authors.
1. The authors should show the viral genome content in the subjects analyzed. 2. I think it is difficult to interpretate the results in this kind of experimental settings since the reference gene expression levels can be affected by several factors not only due to SARS-CoV-2 infections. For example, the age and health conditions of the subjects in the groups can be different. 3. Authors did not mention the quantity of cDNA was used in the RT-PCR reactions. 4. It was mentioned that the groups were categorized by institutional definitions, but the categorization of the groups must be approved by certain international bodies otherwise, it is hard to apply the findings of this study in general. 5. In Fig.2, some of the ct values for the GUSB or PGC-1a or GAPDH in pre-CAM or post-CAM are above ct 30.0, in some cases even above ct 34.0. I think the ct values above 30 is not reliable for calculations. Moreover, the PCR was performed at annealing temperature of 60 degrees for all the reference genes. However, the primers could have better amplification efficiency at other temperatures as well. Not sure, if authors have performed gradient PCR and selected the annealing temperature for the primers. 6. A 260/280 ratio of ~2.0 is generally accepted as "pure" for RNA. In Table 1, the RNA quality for the groups of Severe, pre-CAM and post-CAM is below 2.0 and varies a lot compared to other groups. This is very important factor in the whole experiment as the downstream steps depends on the purity of RNA used for cDNA preparation. 7. To support the RT-qPCR data, authors should have analyzed the expression of these reference genes by Western blotting from the PBMC's collected. 8. Based on normalization with CypA, NRF2 expression level did not correlate with the disease severity in the groups. Therefore, authors should have analyzed for at least 4-5 genes that may have differential expression in SARS-CoV-2 infections to prove that CypA is the suitable reference gene compared to GAPDH.

Dear Dr. Toth,
At the outset, we would like to thank the reviewers for taking out time and providing valuable comments that has allowed us to further strengthen the manuscript. In the light of these comments, we have performed new experiments including gene expression of two additional target genes (IL-6 and IL-15), supporting cytokine (IL-6) and protein quantification data (CRP), and rephrased interpretations in the light of new data that has further strengthened our conclusions. Therefore, we believe that in the revised manuscript, we have adequately addressed the concerns raised by the reviewers. A point-by-point response is attached below: Reviewer #1  (10), 617-621. We have now added statements in the revised manuscript alluding to these facts (Lines 251-53; and Lines 351-55).
Comment 2: I think it is difficult to interpretate the results in this kind of experimental settings since the reference gene expression levels can be affected by several factors not only due to SARS-CoV-2 infections. For example, the age and health conditions of the subjects in the groups can be different.

Response:
We don't understand "this kind of experimental" settings as this work involves patients' samples, and the settings are as good as for any clinical study. The subjects in our data set represent the real-world patient heterogeneity which is not adequately captured by the standard laboratory experimental set-up including with the animal studies. The whole premise of the study was to select an appropriate reference gene that could capture precisely these kind of broad-spectrum manifestations in the patients having COVID-19 infection (asymptomatic, symptomatic, mild, moderate and severe) including co-infection of mucormycosis. The very definition of the ideal reference gene warrants that the gene expression should be minimally regulated across the disease spectrum so as to have minimal influence of the patient heterogeneity in terms of age/sex or severity, treatment etc. It is one of the primary rationales to experimentally determine an appropriate reference gene.
In this context, we screened nine popularly used reference genes as candidates and determined how their expression vary across patient groups. These patient groups were manifesting different clinical symptoms and had a broader age profile (Table 3 original manuscript) even though the mean age of the patients was similar across the groups ( Figure  1, rebuttal).

Figure 1: Mean age of patients enrolled in particular group for this study.
Graph represent Mean±SEM. One-way ANOVA was used to ascertain significance and no significant difference exist among various groups.
These patients received different treatments as per the clinical requirement ascertained by attending physician based on risk and severity (asymptomatic, symptomatic, mild, moderate and severe), and co-infection of mucormycosis. As per the human ethical approval guidelines, we have clearly stated in line no 383 that "The design of study did not influence the clinical care or treatment for any of the patient." Indeed, a large number of the nine candidate reference genes displayed variations in gene expression based on patient heterogeneity as was expected. That is why those genes did not satisfy the prerequisites of a suitable reference gene. Overall, only CypA demonstrated consistent expression across all the patient groups irrespective of any possible clinical confounder and accordingly emerged as the most suitable reference gene in this study. Our work is therefore consistent with the guidelines to reference gene selection for quantitative real-time PCR (Radonić et al, 2004) as the reference genes identified in our study have expression prevalence across the PBMCs of various categories of patients including that of co-infection, and are minimally regulated allowing the accuracy of RNA transcription analysis. The study design also adheres to the established standards of clinical research including the one used to determine reference genes in both humans and animals (Dheda et al, 2004;Montero-Melendez & Perretti, 2014;Pilbrow et al, 2008;Silver et al, 2006).
Based on the scientific arguments presented above, we believe that our study truly captures the real clinical scenario and highlights the heterogeneity of expression in popular/commonly used reference genes; and finally identifies CypA as an ideal reference gene whose expression remains consistent across the clinical spectrum. Wherever necessary we have now expanded on our discussion emphasizing on these facts. We thank the reviewer for the comment as it allowed us to emphasize more clearly on the limitations of conventional reference genes and suitability of CypA as reference gene.

Comment 3: Authors did not mention the quantity of cDNA was used in the RT-PCR reactions.
Response: We apologise for this oversight. The final concentration of cDNA was 10 ng/reaction in qRT-PCR and the same is now included in the revised manuscript (Line no 416). Comment 4: It was mentioned that the groups were categorized by institutional definitions, but the categorization of the groups must be approved by certain international bodies otherwise, it is hard to apply the findings of this study in general. Response: Thank you for highlighting this important point. AIIMS, New Delhi (Our institute) was the National Nodal Centre for COVID-19 and follows the same guidelines as recommended by Indian Council of Medical Research, ICMR, and Govt. of India.
These guidelines are also similar to guidelines from other international organisations such as WHO, and NIH, USA. In the revised work, we have appropriately rephrased and expanded on these guidelines (lines no 363-83). Fig.2, some of the ct values for the GUSB or PGC-1a or GAPDH in pre-CAM or post-CAM are above ct 30.0, in some cases even above ct 34.0. I think the ct values above 30 is not reliable for calculations. Moreover, the PCR was performed at annealing temperature of 60 degrees for all the reference genes. However, the primers could have better amplification efficiency at other temperatures as well. Not sure, if authors have performed gradient PCR and selected the annealing temperature for the primers. Response: Thank you. We respond to this comment in two parts as indicated below: a) Some of the ct values for the GUSB or PGC-1a or GAPDH in pre-CAM or post-CAM are above ct 30.0, in some cases even above ct 34.0. I think the ct values above 30 is not reliable for calculations.

Comment 5: In
We appreciate the reviewer's comment. First, we respectfully disagree with the subjective opinion that values above Ct 30 are not reliable. This assertion belies scientific evidence, and contradicts the well-established and universally accepted scientific premises/principals. Essentially, Ct values represent transcriptional level of a gene in that particular cell/condition; if any gene is less represented/less abundant, that will be reflected in the Ct values. Based on the natural expression, the level may be less or high which will be reflected in the Ct value. Ct values actually are considered in acceptable range until it does not cross the limit of detection or negative control/NTC (Non-template control) and rather than arbitrary cut-off values as indicated in the comment. Guidelines to reference gene selection for quantitative real-time PCR recommends a Ct ≤40 for determination of qRT-PCR efficiency and intra-and inter-assay variability for selection of appropriate gene (Radonić et al., 2004). Consequently, literature is replete with the studies supporting the Ct values above 30 for scientific interpretations and experimentations including confirming of clinical status and follow-up care (Dheda et al., 2004;Radonić et al., 2004;Silver et al., 2006). Heterogeneity in gene expression among patients is expected, and commonly present, and same is aptly reflected in the Ct values. It is precisely this heterogeneity that makes it essential to determine the most suitable reference gene, especially when working with patients. We have adequately highlighted these facts by providing ranges of Ct values observed in our study (Lines 108, lines 119-124, original manuscript). In our manuscript, the differences in Ct values of a specific reference gene from a particular group that the reviewer is alluding to are in fact natural variations in expression observed among patients (Brym et al, 2013;Gao et al, 2008;Silver et al., 2006;Yu et al, 2008). It is well reported that these housekeeping genes are involved in plethora of cellular pathways, which may get dysregulated in patients based on variety of factors including severity (Dheda et al., 2004;Montero-Melendez & Perretti, 2014;Shivashankar et al, 2021). Further, most qRT-PCR based COVID-19 diagnostics kits approved by WHO, CDC and other National and international agencies consider acceptable range is ~ Ct ≤40 (Table 1, rebuttal). In our study, irrespective of variations in severity of infection, the overall Ct values remain within the 5 range of each group (dispersal range 6.43-39.69; lines 108 original manuscript) and we could successfully identify the most stable reference gene to be used during COVID-19. Thank you. Firstly, Figure 1 of our original manuscript clearly show a single melt curve and melt peak for each individual gene. This is indicative of a single amplicon denoting specificity of the primers at the conditions used. Further, primer designing was performed in a manner to have their expected best amplification efficiency at 60°C; and it was validated using healthy subjects. Nonetheless, to address reviewer's concern, we again checked amplification efficiency of primers of individual genes through gradient PCR by selecting a temperature range of 58°C-61°C and with three different subjects (Figure 2, rebuttal).
As shown below, we observed the best amplification efficiency for primers for each candidate reference gene at 60°C only. Comment 6: A 260/280 ratio of ~2.0 is generally accepted as "pure" for RNA. In Table  1, the RNA quality for the groups of Severe, pre-CAM and post-CAM is below 2.0 and varies a lot compared to other groups. This is very important factor in the whole experiment as the downstream steps depends on the purity of RNA used for cDNA preparation.

Response:
We understand the reviewer's concern. 260/280 ratios actually can be impacted by several factors including the complexity of clinical samples, storage, template composition, and even solution used to suspend RNA (Okamoto & Okabe, 2000;Willinger et al, 2017;Yu et al., 2008). While 260/280 ratio of ~2.0 is desirable, the 260/280 values ~1.7 are acceptable and frequently used (Ghawana et al, 2011;Willinger et al., 2017;Zhang & Yang, 2022 In our study, comparison of 260/280 ratio across all groups revealed that there were no significant differences in the mean 260/280 ratio among various groups including the severe COVID-19 and CAM groups (Figure 3, rebuttal). There was no significant difference in the 260/280 ratio (Y-axis) across the groups (X-axis). Data represent mean±SD; One-Way ANOVA with Kruskal-Wallis test with post-hoc corrections. "ns" represents nonsignificant differences (*, P-value ≤0.05). This figure is now added as supplemental Fig. S2 in revised manuscript.
The mean 260/280 ratio was in the range of 1.94 to 1.98) for healthy, asymptomatic and mild infected; and (1.73 to 1.79) in case of CAM and severe COVID-19 patients. Incidentally, the mortality rates among these cohorts of severe COVID-19 and CAM were over 30%; and ~50%, respectively. Nonetheless, to rule out the effect of these-minor and statistically nonsignificant variations in 260/280 ratio among various groups, we re-analyzed our data by excluding patient data with 260/280 ratio <1.7, we observe that this did not alter the conclusions of our analysis and CypA remained the most suitable reference gene.
In the revised manuscript, we have added points pertaining to this discussion (Lines no 292-99) Comment 7: To support the RT-qPCR data, authors should have analysed the expression of these reference genes by Western blotting from the PBMC's collected.

Response:
We respectfully disagree. The aim of our study was to find the most suitable reference gene for qRT-PCR data analysis in COVID-19 and CAM patients across the disease spectrum. It is well known that qRT-PCR data measures expression of gene (transcription), while the WB analysis incorporates additional elements including mRNA half-life, translational efficiency, trafficking, stability and post-translational modifications.
Furthermore, SARS-CoV-2 is known to affect splicing, translation and trafficking of host proteins to suppress host response (Banerjee et al, 2020;Finkel et al, 2021). During infection, SARS-CoV-2 develops multiple strategies like degradation of host mRNA pool to hijack host translation machinery to synthesize its own protein (Finkel et al., 2021;Fisher et al, 2022). Considering these well-established premises, it is not clear to us how WB analysis would have supported meaningful interpretation of qRT-PCR data. This fact is also clearly reflected in the similar studies pertaining to identification of suitable reference gene for qRT-PCR data analysis in other organisms/conditions (Table 2, rebuttal), none of which have performed WB analysis.  Nonetheless, in the revised manuscript, we have now determined the IL-6 gene expression using normalization with CypA, and GAPDH. We correlated gene expression profile with the levels of cytokine IL-6 in these patients. We demonstrate that gene expression values for IL-6 obtained by normalization with CypA, and not with GAPDH, are consistent with the IL-6 levels in plasma. These new data and discussion are now added in the revised manuscript (Lines no 203-216; and 269-71). Comment 8: Based on normalization with CypA, NRF2 expression level did not correlate with the disease severity in the groups. Therefore, authors should have analysed for at least 4-5 genes that may have differential expression in SARS-CoV-2 infections to prove that CypA is the suitable reference gene compared to GAPDH.

Response:
We thank the reviewer for this insightful comment. To confirm the findings of CypA based normalization vis a vis GAPDH based normalization, we wanted to use the genes whose expression is previously confirmed in COVID-19 patients using RNA sequencing. This would allow us to evaluate how closely qRT-PCR data normalized using CypA or GAPDH compares with the RNA sequencing-based studies. It is in this context that we used major redox regulator NRF2 whose expression is dysregulated during COVID-19, and is not sufficiently up-regulated, and remained relatively suppressed/minimally induced to mount an effective anti-oxidant response . During oxidative stress conditions such as inflammation, NRF2 is usually significantly induced to maintain redox homeostasis. The relative suppression of NRF2 expression by SRAS-CoV-2 is a viral strategy to hijack host protective responses (Olagnier et al, 2020; and induction of NRF2 pharmacologically is proposed to be an effective strategy against COVID-19 (McCord et al, 2021) .
We first show that our data on NRF2 expression based on CypA-normalization is actually consistent with the previous reports on NRF2 expression showing that SARS-CoV-2 infection dysregulates NRF2 to impair its proper induction leading to an increased oxidative stress, a hallmark feature of COVID-19 (Olagnier et al., 2020;. Similar to RNA sequencing-based study, normalization with CypA revealed a relatively low NRF2 expression even in the oxidative stress environment of COVID-19-associated inflammation. Moreover, we did not observe significant difference in the expression of NRF2 across the spectrum of COVID-19 infection (Figure 4, rebuttal and Fig. 4 in revised manuscript), as is also indicated from RNA Sequencing analysis (Jain et al., 2021). On the other hand, normalization with GAPDH in our study rather indicated massive up-regulation of NRF2 expression except in case of severe COVID-19 patients. NRF2 expression based on CypA-normalization reveals no significant difference exists in its expression across various patients' groups. This is consistent with the RNA sequencing-based expression wherein NRF2 expression is relatively low even in the high oxidative stress environment reported for COVID-19 and CAM infections (Jain et al., 2021;Olagnier et al., 2020). On the other hand, normalization with GAPDH reveals significant expression heterogeneity that varies from 2-fold to over 600-fold across spectrum of COVID-19 and CAM infections. This does not match with RNA sequencing-based data from published works and erroneously creates the impression of NRF2 over-expression. Healthy subjects, represents the baseline NRF2 expression. Data represents mean ± standard error of the mean (SEM) performed with at least eight subjects in duplicates. Significance is ascertained by using a One-way ANOVA with Kruskal-Wallis test with Dunn's post-hoc corrections, N=8/group (*, P-value <0.05; **, P-value < 0.01; ***, P-value < 0.001, and ns=non-significant).
For further strengthening of our conclusions about utility of CypA, in the revised work, we have additionally determined the expression of IL-6 and IL-15 genes ( Figure 5 and 6 of rebuttal, and Fig. 5-7 in revised manuscript). Furthermore, the data on IL-6 gene expression, a marker for inflammation is also corelated with cytokine IL-6 levels. The inflammation status as suggested by IL-6 is further corroborated by determining CRP levels.
IL-6 gene expression: IL-6 expression is regulated by NRF2 as NRF2 blocks the proinflammatory responses by inhibiting the expression of IL-6 at the transcription level (Kobayashi et al, 2016;Matsuoka et al, 2016). A negative relationship therefore exists between expression of NRF2 and IL-6, and a high NRF2 expression should lead to suppression of IL-6 expression. Surprisingly, we observe that the normalization with GAPDH as reference gene produces a high NRF2-high IL-6 expression phenotype as both IL-6 and NRF2 showed massive induction (>400 fold) in case of patients with mild, and CAM (pre and post) infections. This shows that GAPDH is not an ideal reference gene for qRT-PCRbased normalization during COVID-19. Moreover, we also did not observe any significant differences in the expression of IL-6 between asymptomatic, and severely infected COVID-19 patients when data is normalized with GAPDH. This is again in contrast to the established fact that inflammation is minimal in asymptomatic, and very high in patients having severe COVID-19. With GAPDH, we observed similar expression among asymptomatic, mild, moderate and severe infections. On the other hand, RNA sequencing data indicates a severity dependent increase in IL-6 expression going from mild to severe infections (Mukhopadhyay et al., 2021;Islam et al., 2021;Delgado-Roche et al., 2020). Only IL-6 expression derived from normalization with CypA, therefore, is in agreement with previously published studies. Y-axis represents differences in the fold-change based on data normalization vis à vis CypA or GAPDH in various groups (X-axis) of COVID-19 and CAM patients. Significance is ascertained with respect to healthy subjects by using a One-way ANOVA with Kruskal-Wallis test with Dunn's post-hoc corrections (*, Pvalue <0.05; **, P-value < 0.01; ***, P-value < 0.001, and ns=non-significant).
On the other hand, normalization of data with CypA showed an overall higher IL-6 and low NRF2 expression in all groups consistent with the established relationship of IL-6 and NRF2 expression. The highest IL-6 expression was observed in severe COVID-19 patients that had the lowest NRF2 expression. Furthermore, expression of IL-6 in asymptomatic, mild and severe cases differed significantly, as it should be as inflammation levels vary significantly among these groups with asymptomatic having the lowest and severe having the highest inflammation.
To further support our observations, we next determined IL-6 levels in plasma (Figure 6, rebuttal, and Fig. 6a in revised manuscript) and noticed that compare to GAPDH, the CypA normalized IL-6 expression concurs well with the IL-6 levels further validating CypA as a reference gene. Asymptomatic patients have IL-6 levels comparable to healthy individuals. However, GAPDH normalized data showed non-significant differences in IL-6 expression between severe and asymptomatic. level. These data support the gene expression profile derived from normalization with CypA. Significance is ascertained by using a One-way ANOVA with Kruskal-Wallis test with post-hoc corrections. N=8-12/group (*, P-value <0.05; **, P-value < 0.01; ***, and P-value < 0.001).

Estimation of CRP level:
To additionally strengthen our conclusions, we determined plasma CRP (C-reactive protein) levels to reflect on the status of systemic inflammation among these patients. We observed a direct correlation between CRP and disease severity, with asymptomatic subjects having the lowest, and severe COVID-19 patients having the highest CRP levels (Figure 7 in rebuttal, and Fig. 6b in revised manuscript). We observed that this data concurs with the plasma IL-6 levels, and IL-6 gene expression when CypA is used for normalization. Significance is ascertained by using a One-Way ANOVA with Kruskal-Wallis test with post-hoc corrections. N=8-12/group (*, P-value <0.05; **, P-value < 0.01; ***, and P-value < 0.001).

IL-15 gene expression:
In addition, we also scored for the gene expression of IL-15 and observe that the expression profile obtained following normalization with CypA corresponds better with the RNA sequencing based published work. Normalization with GAPDH does not concur with RNA-sequencing-based data which shows that IL-15 expression increases significantly from mild to moderate, and then returns to baseline expression in severely sick (Jain et al., 2021). In case of normalization with CypA, expression pattern concurs with RNA sequencing data by demonstrating an increase from mild to moderate, and a reduction in severely sick. Y-axis represents differences in the fold-change based on data normalization in various groups (X-axis) of COVID-19 and CAM. Significance is ascertained by using a One-way ANOVA with Kruskal-Wallis test with Dunn's posthoc corrections (*, P-value <0.05; **, P-value < 0.01; ***, and P-value < 0.001).
Taken together, our additional data involving expression profiling of genes namely IL-6, and IL-15 along with determination of plasma levels of IL-6 and CRP across the spectrum of COVID-19 and mucormycosis infections firmly establish the utility of CypA as the reference gene.
The statements alluding to the above facts and additional data have now been added to the revised manuscript (lines no 190-221 and 302-336).

Reviewer #2 (Comments for the Author)
The

Strength
Studying the co-infection of SARS-CoV-2 and mucormycosis is relevant based on existing public health issue and the lack of studies published in the field. In addition, the selection of accurate reference genes in different disease outcomes could advance the field.

Weaknesses
The observations postulated in the paper are novel and interesting. However, the discussion is not robust enough to make convincing conclusion. The authors should show more data from published studies related to their results.

Response:
We thank the reviewer for the encouraging words and highlighting strength and weaknesses.
In the revised manuscript, we have now significantly strengthened our conclusions by providing additional data, including more candidate genes for validation of CypA as reference gene and supplementing this data by quantifying levels of CRP and IL-6 in the plasma of the patients. In the light of additional data, and a significantly improved discussion, we believe that we have now adequately addressed the comments in the revised manuscript. We thank the reviewers for their time and efforts to improve our work. Response: We appreciate the reviewer for this comment and apologise for this confusion, which has arisen primarily due to Figure 3b. Following the reviewer's comment, we revisited this figure and realised an inadvertent error in the processing of Bestkeeper data that factored in only coefficient of correlation (r 2 ) without incorporation of the stability values. Now, we have fixed this oversight and reprocessed this data (Fig. 3b in revised manuscript, and Figure 9 of rebuttal), where CypA showed highest stability value whereas GAPDH showed least stability value. Therefore, the stability value of CypA is now in accordance to Figure 2 of original manuscript. We thank the reviewer for bringing this to our notice that has allowed us to fix it. (a) Intra-group average standard deviations for individual patient groups were obtained using delta Ct method, and are added to obtain a cumulative average standard deviation (Y-axis). The highest cumulative average standard deviation (i.e., GAPDH) indicates the least gene stability (X-axis). (b) The graph represents a cumulative stability value (Y-axis) obtained by adding intra-group stability value based on coefficient of correlation (r) and standard deviation using Bestkeeper algorithm. The candidate gene yielding the lowest cumulative value (i.e., PGC-1α) indicates the highest stability. (c) In Normfinder algorithm, intra-group stability value (Y-axis) based on coefficient of variations are added to obtain a cumulative stability value. The candidate gene yielding lowest cumulative stability value (i.e., CypA) indicates the most stable reference gene. (d) In GeNorm algorithm, intra-group expression stability (M value) based on the average pairwise variation of an individual gene with all other genes used as control to obtain a cumulative expression stability (Y-axis). The gene yielding the highest M value (i.e., GAPDH) indicates the lowest stability. (e) In RefFinder algorithm, intragroup geometric mean (Y-axis) based on geometric mean of ranking values calculated by all above algorithms are added to obtain a cumulative geometric mean of ranking values. CypA yielded the lowest cumulative geometric mean indicating the highest stability and therefore emerged as the most suitable reference gene.

Major comments
Furthermore, GAPDH is a multi-faceted gene, which in addition to glycolysis, is also involved in diverse cellular functions including regulation of inflammation (Montero-Melendez & Perretti, 2014). Inter-group variability in GAPDH expression can also be due to the differential inflammation status observed in COVID-19 patients as is also highlighted by IL-6 and CRP data presented in the revised manuscript (Millet et al, 2016;Shivashankar et al., 2021). We have now expanded these statements in the discussion section of the revised manuscript (Lines no-268-271).
Comment 2: Figure 4: As the health subjects represent the baseline NRF2 expression, is there significant difference in NRF2 expression among samples from moderate, asymptomatic and/or pre/post-CAM diseases when CypA was used as internal control gene? If so, how can the NRF2 expression variability among them using CypA as internal control gene be explained?

Response:
We appreciate the reviewer's comment. In the revised manuscript, we have performed the statistical analysis with the data of NRF2 gene expression from healthy subjects as baseline. Even though there are differences in expression of NRF2 across various categories of patients, these differences are not statistically significant when data is normalized with CypA as internal control. This is consistent with the earlier reports in COVID-19 where in NRF2 expression is relatively suppressed/minimally induced even though the inflammation levels and consequent oxidative stress levels are very high (Olagnier et al., 2020;. NRF2 is a major regulator of host anti-oxidant response, which also regulates the transcription of cytoprotective genes to mitigate the oxidative stress mediated tissue damage, and inflammatory response induced by viral infection. Nonetheless, during COVID-19, NRF2 expression is dysregulated and not sufficiently induced to mount an effective anti-oxidant response . The relative suppression of NRF2 expression by SRAS-CoV-2 is a viral strategy to hijack host protective responses (Olagnier et al., 2020; and induction of NRF2 expression is proposed to be an effective strategy against COVID-19 (McCord et al., 2021). The slightly higher NRF2 levels observed in CAM group can be a response to mucormycosis infection, and NRF2 expression returns closer to base-line following reduction in fungal burden post-surgery (post-CAM).
We have alluded to these facts in the revised manuscript. Also, please refer to our response to comment 8, reviewer 1.
Comment 3: Also, at line 164, the authors comment that NRF2 is suppressed in Covid-19 patients, which contradicts the Figure 4 data that shows NRF2 expression seems to be higher mainly in Moderate COVID-19 and pre-CAM disease comparing to health subjects even using the CypA as the internal control. Thus, the authors should specify the difference in the expression profile of NRF2 inter-group based on the usage of CypA.
Response: As explained above, we have used NRF2 only as a candidate gene whose expression was previously reported in COVID-19 patients using RNA sequencing. We now provide data using additional candidate genes namely IL-6 and IL-15 to support CypA as internal control gene. As far as NRF2 is concerned, multiple studies have shown that expression of NRF2 is suppressed/ minimally induced in COVID-19 and that regulation of NRF2 expression by SARS-CoV-2 is a viral strategy to hijack host protective responses . We first show that our data on NRF2 expression based on CypA-normalization is actually consistent with the previous reports on NRF2 expression showing that SARS-CoV-2 infection dysregulates NRF2 to impair its proper induction leading to an increased oxidative stress, a hallmark feature of COVID-19 (Olagnier et al., 2020;. Similar to RNA sequencing-based studies, normalization with CypA revealed a relatively low NRF2 expression even in the oxidative stress environment of COVID-19-associated inflammation. Moreover, we did not observe significant difference in the expression of NRF2 across the spectrum of COVID-19 infection (Figure 4, rebuttal and Fig. 4 in revised manuscript), as is also indicated from RNA Sequencing analysis (Jain et al., 2021).
There is no published study involving CAM patients for us to correlate or compare NRF2 expression. Slightly higher NRF2 levels observed in CAM group in our study can be a response to mucormycosis infection, and NRF2 levels returns closer to base-line following reduction in fungal burden post-surgery (post-CAM) (Figure 4 of rebuttal, and Fig. 4 in revised manuscript). Nonetheless, further focused studies would be required to study NRF2 and redox regulation in CAM patients which is beyond the scope of this work.
We have incorporated these comments in the revised manuscript.  (Table 2). Response: We apologise for this confusion and this have been fixed in the revised manuscript.

Minor comments
Taken together, in the revised work, we have added significant amount of new data including gene expression and protein estimation to adequately justify the relevance of CypA as the most suitable gene for qRT-PCR based gene expression analysis.
We thank both the reviewers for valuable comments and constructive suggestions that has allowed us to improve the manuscript significantly.